##genesurvival


#引用包
library(survival)
library(survminer)

inputFile="rs_tcga.txt"      #风险文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\13.genesurvival")     #设置工作目录

##提取模型基因表达量
rt=read.table("OV_TPM.txt", header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#读取基因列表文件，提取免疫检查点相关基因的表达量
gene=read.table("modelgene.txt", header=T, sep="\t", check.names=F)
sameGene=intersect(row.names(data), as.vector(gene[,1]))
data=t(data[sameGene,])
data=log2(data+1)

#删除正常样品
group=sapply(strsplit(row.names(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[group==0,]
row.names(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(data))
data=avereps(data)
#读取输入文件

rt=read.table(inputFile, header=T, sep="\t", check.names=F, row.names=1)
rt=rt[,1:(ncol(rt)-1)]
rt=merge(rt,data,by = 0)
rownames(rt)=rt$Row.names
rt=rt[,-1]

#对基因进行循环
outTab=data.frame()
for(gene in colnames(rt)[3:ncol(rt)]){
  if(sd(rt[,gene])<0.001){next}
  data=rt[,c("futime", "fustat", gene)]
  colnames(data)=c("futime", "fustat", "gene")
  
  #获取最优cutoff
  res.cut=surv_cutpoint(data, time = "futime", event = "fustat", variables =c("gene"))
  res.cat=surv_categorize(res.cut)
  fit=survfit(Surv(futime, fustat) ~gene, data = res.cat)
  
  #比较高低表达生存差异
  diff=survdiff(Surv(futime, fustat) ~gene,data =res.cat)
  pValue=1-pchisq(diff$chisq, df=1)
  outVector=cbind(gene, res.cut$cutpoint[1], pValue)
  outTab=rbind(outTab,outVector)
  if(pValue<0.001){
    pValue="p<0.001"
  }else{
    pValue=paste0("p=",sprintf("%.03f",pValue))
  }
  
  #绘制生存曲线
  surPlot=ggsurvplot(fit,
                     data=res.cat,
                     pval=pValue,
                     pval.size=6,
                     pval.method = TRUE,
                     legend.title=gene,
                     legend.labs=c("high","low"),
                     xlab="Time(years)",
                     break.time.by=1,
                     palette = c("#CB4535","#273359"),
                     ggtheme = theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+
                       
                       theme(panel.grid = element_blank()),
                     conf.int=F,
                     title = "Survival curve",
                     risk.table=F,
                     risk.table.title="",
                     risk.table.height=.25,
                     font.main = c(16, "bold.italic", "black"), # 设置标题字体大小、格式和颜色
                     font.x = c(14, "bold.italic", "black"), # 设置x轴字体大小、格式和颜色
                     font.y = c(14, "bold.italic", "black"), # 设置y轴字体大小、格式和颜色
                     font.tickslab = c(12, "plain", "black"),# 设置坐标轴刻度字体大小、格式和颜色
                     linetype = "solid",#改变线条的渐变
                     legend=c(0.8,0.8))##图例位置
  
  #输出图形
  pdf(file=paste0("Survival.",gene,".pdf"), width=4.5, height=4, onefile=FALSE)
  print(surPlot)
  dev.off()
}
